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We perform numerical simulations of the critical gravitational collapse of a massive vector field. 
The result is that there are two critical solutions. One is equivalent to the Choptuik critical solution 
s ! ' for a massless scalar field. The other is periodic. 

>V 



co 



in 
o 

CO 



5h . 

tofj 



5-H 



I. INTRODUCTION 



Critical gravitational collapse was first found by Choptuik^J in simulations of a spherically symmetric massless 
I | scalar field. A natural question to pose is then how critical collapse behaves when the scalar field has a mass, since 
. this will introduce a characteristic length that destroys the scale invariance of the field equations. This question 
t^J- 1 was studied by Brady et. aZ.Q The results of reference^] show that there are two critical solutions: one which is 
essentially the Choptuik critical solution for the massless scalar field, and another which is a periodic solution first 
found by Seidel and Suen . At first it might seem puzzling that the Choptuik solution can be a critical solution for 
both the massless and massive scalar field. The resolution of this conundrum is that as the singularity is approached 
in the Choptuik critical solution, the amplitude of the scalar field remains bounded while its gradient diverges. In the 
stress energy tensor, the mass terms are associated with the amplitude of the field, while other terms are associated 
■ — . with its gradient. So as the singularity is approached the mass terms in the stress energy become negligible. 

Given the results of reference^ one might conjecture that similar behavior occurs in the case of a spherically 
symmetric massive vector field: i. e. that there is a critical solution for which the mass of the vector field becomes 
negligible. However, this conjecture involves a paradox: a massless vector field is just a Maxwell field, and a spherically 
symmetric Maxwell field has no degrees of freedom. Therefore there is no gravitational collapse (and thus no critical 
solution) of a spherically symmetric massless vector field. What then is the critical behavior of a massive vector field? 



In this paper we consider this question. We perform numerical simulations of the collapse of a spherically symmetric 
massive vector field. The equations and numerical methods are presented in section 2. Results are given in section 3 
and conclusions in section 4. 

II. EQUATIONS AND NUMERICAL METHODS 

A massive vector field is described by the Proca Lagrangian 

C=-^F ab F ab -^ 2 A a A a (1) 

where F a b = V a Ai, — VbA a and fi is a constant. Note: throughout, we use the conventions of Wald£| and in particular 
use a metric with signature (— , +, +, +). Note that were we to employ the opposite signature, the sign of one of the 
terms in the Proca Lagrangian would have to be changed (see e.g. 5]). The equation of motion that follows from 
equation Q is 

V a F ab - fi 2 A b = (2) 

from which it follows that 

V a A a = (3) 



2 



It also follows from equation (JTJ that the Einstein field equation is 



Gab = 2F ac F b c + 2n 2 A a A b - g ab 



2 FcdFC 



(4) 



We now specialize to spherical symmetry. We employ two different methods to simulate the Einstein-Proca sys- 
tem: a Cauchy method using polar-radial coordinates and a characteristic method using the coordinates used by 
Christodouloujg to treat the Einstein-scalar system. The metric in polar-radial coordinates takes the form 



ds 2 



-a 2 dt 2 + a 2 dr 2 + r 2 (d0 2 +: 



Define the quantities X and W by 



X 
W 



-A t 



— (d t A r 
aa 
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From the definition of W we find 



Then from equation @ we have 



d t A r = aaW - 
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From equation J2J) we find the following constraint equation for W 

\d r {r 2 W) +fi 2 X = 



(5) 

(6) 
(7) 

(8) 
(9) 
(10) 



Equations (JHJ), © and QlOjl provide the evolution equations for the matter field. In spherical symmetry the metric 
has no degrees of freedom. Therefore the metric functions a and a are given by "constraint" equations once the 
matter fields are known. To find the appropriate constraint equations, note that for a metric of the form given in 
equation JSJ) the tt and rr components of the Einstein tensor are 



G tt = - 
i 

It then follows from equations (|llf> and Q that 

^ 2 1 



, n da 
a + a A + 2r— 
or 

, „ da 

1 - a 2 )a + 2r— 
or 



Ida 
a dr 



^—±= r -[a 2 W 2 + n 2 (X 2 + A 2 ,)] 



Then using equations Qllfl. (|12|l and 10} we find 



1 da I da , „ 
a dr a or 
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(11) 
(12) 

(13) 
(14) 



To implement these equations numerically, we replace spatial derivatives with centered differences and implement the 
time evolution using the iterated Crank-Nicholson method. We also put in Kreiss-Oliger dissipation 8] for added 
stability. Initial data for this system is X and A r at the initial time. Given these initial data, equations (|l(jp. \YA\ 
and (|14|) are then integrated in turn to obtain W and the metric functions. Finally equations J3J) and are used to 
produce A r and X at the next time step. 

In addition to this unigrid code, we also perform simulations with an adaptive mesh code. Note that our equations 
are quite similar to those used by Choptuik et al |9( to simulate critical collapse in the Einstein- Yang-Mills system. 
Our adaptive code is produced by modifying the code of reference to simulate our system. 

We now present the characteristic method using the coordinates of reference [fj. Here the metric takes the form 



ds 2 



e 2v du 2 - 2e v+A dudr + r 2 dn 2 



(15) 
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We introduce the null vectors 



I" 



d_ 

On 



±e" A 



dr) 

d_V 

dr . 



(16) 
(17) 



The matter in this coordinate system is determined by the components A u , A r and F ur — e v+ W. Note that W 
defined in this way is the same as in equation(|7J) as can be seen by the fact that F ab F a i, — —2W 2 in both cases. We 
also introduce the quantities g = e u+x and g = e v ~ x . 
For a metric of the form of equation i|15|) we have 



G ab n a l b = 

gr z 



n ,a,b _ 2 9 dg 
rg z or 

or 



From equation 10} it then follows that the corresponding Einstein equations become 

rg z or g 



gr 



= w 2 



(18) 
(19) 

(20) 
(21) 



We now define a scalar field <j> by d r cj> = \iA T . Note that this defines cj> up to addition of an overall constant, since 
smoothness of 4> implies that d r (j> = d u (p at r = 0. We also define the quantities h and h by 



h 
h 



d_ 

dr 
1 



Then the solution of equations (|20|l and l|21ll become 



exp 

1 

r 



hdr 



-{h-hfdr 
r 



g(l - W 2 )dr 



Next we find an expression for the matter variable W in terms of h. Contracting equation (0) with l a we find 

d 



dr 



(r 2 W) + ^r 2 A r = 



for which the solution is 



W=-4/ r(h-h)dr 



We now find an evolution equation for h. From F a i, = 2di a A^ it follows that 

d u A r - d r A u = gW 

Which from the definition of <f> leads to 

d u d r (t) = d r (fiA u ) + figW 
Now from this equation, its integral with respect to r and equation (|22() we find 



d u h = d r (r[iA u ) + fj,grW + / [igWdr 



(22) 
(23) 

(24) 
(25) 

(26) 

(27) 

(28) 
(29) 

(30) 
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Thus to find the evolution equation for h we must find an expression for d r (r[iA u ) in terms of h. To do this, we note 
that from equation it follows that 

d u A r + ^d r [r 2 {A u - gA r )) = (31) 

Then subtracting equation 128|) from equation 1311) we obtain 

d r (nrA u ) = -^fxrgW + ^-d r (r 2 gd r (b) (32) 

Note that equation i|32|) can be integrated to yield an expression for the remaining matter variable A u in terms of h. 
Thus, given h at a time u, all matter and metric variables at that time can be determined by integrals. Now using 
equation and equation in equation (|3l7|l we obtain 

Dh=^(h-h)(g-g)- \grW 2 + \^grW + f figWdr (33) 
2r 2 2 J 

Here D = d u — (g/2)d r is derivative along the ingoing null direction. 

The numerical method used is the same as that used in reference 0] for scalar field collapse. Initial data is given 
for h at an initial u. Equations l|23|) . I|27|) . (|24|l and l|25|l are then integrated to find the other matter and metric 
variables. The integration method is Simpson's rule for unequally spaced points, but near the origin a Taylor series is 
used. Then equation is used to find h at the next value of u. Each grid point is an ingoing light ray, and both h 
and r are evolved along the grid points. When a grid point reaches r = 0, it is removed from the computational grid 
and when half of the points have been lost, they are put back in between the remaining points by using interpolation. 

The scale invariance (A a , g a \,) — » (kA a , k 2 g a b) of the Einstein-Maxwell system extends to the Einstein-Proca system. 
Specifically, if (A a ,g a b,fi) is a solution of equations J2J and |@J then (kA a ,k 2 g a b, n/k) is also a solution of these 
equations, where k is any positive constant. This allows us to set n = 1 without loss of generality, which we do in 
all runs. Note that large k can render ijl negligible while retaining the effect of the additional Proca terms. We will 
encounter this effect in our investigation of the critical behavior of this system. 



III. RESULTS 



Runs were done on various Unix and Linux workstations and on PCs. To test universality, we tried initial data of 
several forms, including a gaussian shape and a sech shape. However, here we will present only the results of runs 
with the gaussian initial data. Specifically, for the Cauchy codes we use the following form of initial data: X = 
and A r = pr exp[— (r — ro) 2 /a 2 ] where p, rg and a are constants. For the characteristic code, the initial data is 
(f> = pr 2 exp[— (r — ro) / 1 a 2 \. In simulations, we fix rg and a and have p as the parameter that is varied. The critical 
value of p (denoted p*) is found by a binary search. 

We find two different critical solutions depending on the value of a. One is a periodic type I critical solution. Figure 
1, produced using the unigrid Cauchy code, shows X at r = as a function of t for this solution. For this run, we 
have cr = 1.5, r = 3.0 and p* = 0.104135195147191. 

The other solution is a type II DSS critical solution which appears to be identical to the Choptuik critical solution for 
a massless scalar field. Figure 2 shows a plot of InM vs \n(p — p*) for solutions above but near the critical one. (Here 
M is the mass of the black hole). This plot was produced using the adaptive Cau chy code. Here, r — 3.0, a = 0.5 
and p* — 0.134075353579. For the Choptuik critical solution, the results of [lTl Il2| show that the graph of In ill 
vs ln(p — p*) is a straight line with a periodic wiggle. Here the slope of the line is called 7 and the period of the 
wiggle is T w = A/(2~f) where A is the period of the DSS critical solution. The simulations of ^2 & ve 7 = 0.374 and 
A = 3.4453 which yields T w — 4.61. We fit the data of figure © to a straight line plus a sine wave. We find that the 
slope of the line is 7 = 0.379 and the period of the wiggle is T w = 4.63. Thus, to the accuracy of our simulation, our 
DSS critical solution gives values of 7 and T w in agreement with those of the Choptuik critical solution. 

Similarly figure 3 shows a plot of In i? max vs ln(p* —p) for solutions below but near the critical one. Here i? max is the 
maximum value of the scalar curvature at the center. This figure was produced using the characteristic code. Here, 
r = 2.0, a = 0.5 and p* = 0.0501805022078927. (Note though that due to the different type of data, these parameters 
have different meaning than in the Cauchy case). As shown in ^| this sort of plot should also be a straight line with 
a periodic wiggle. Here the slope of the line should be —27 and the period of the wiggle should be T w — A/ (27) A 
fit of the data of this figure to a straight line plus a sine wave yield that the slope of the line is —27 = —0.727 which 
gives rise to 7 = 0.363 while the period of the wiggle is T w — 4.64. These values are again comparable to the values 
of the Choptuik critical solution. 
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FIG. 1: X(0) for the periodic critical solution 




FIG. 2: InM us ln(p — p*) near the DSS critical solution 



We now make a direct comparison between this critical solution and the Choptuik critical solution for a massless 
scalar field. The simplest way to do this is to note that our characteristic equations I|22l - I25(l . (|27[) and i-VMi formally 
go over to the corresponding equations for a massless scalar field if we set the parameter /z to zero. We will return to 
this point in the next section. Thus we can find the Choptuik critical critical solution with our code by performing a 
binary search with = 0. Figure 4 contains a comparison of the two critical solutions. What is plotted is h at r = 
as a function of T where e~ T = u* — u and u* is the value of u at which the singularity forms. We use the invariances 
of the two systems to choose offsets in T and h so that the two solutions coincide at a maximum of h. Note that the 
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two solutions (after an initial transient has died away) are the same. 

IV. CONCLUSIONS 



Given the results of reference 2] for a massive scalar field, it is not surprising that a massive vector field has a type 
I periodic critical solution. What does seem surprising is that it has a DSS critical solution that is identical to that of 
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a massless scalar field. For a DSS critical solution we would expect that since length scales are becoming arbitrarily 
small, that /x is becoming negligible compared to the inverse of the relevant length scale (or rather, since we set [i = 1, 
that the relevant inverse dimensionless length scale is becoming arbitrarily large). 

Thus the DSS critical solution should in some sense also be a solution of the "/i — > limit" of the equations. We 
have already seen how to make sense of this limit in the case where spherical symmetry is imposed and the system is 
expressed in terms of variables chosen to be similar to those of reference . We now show how to make sense of this 
limit more generally using equations (0-EJl. If A a itself has a smooth fi — > limit, then the ^ — ► limit of equations 
(J2I-BJ is simply the Einstein-Maxwell equations. Instead we assume that A a takes the form 

A a = -Pa + HQa (34) 

A* 

Then in order that the stress-energy have a non-singular fi — > limit, we must have V[ a f^i =0 and therefore there 
must be a scalar field (j) such that P a — V a </>. Then in the fi — ► limit equations Q and @ become respectively 

V Q V a = (35) 
Gab = 2V a 0V 6 - 5afc V c 0V c (36) 

In other words, the [i — > limit of the Einstein-Proca system becomes the Einstein-scalar system. 

Consequently it is not surprising that these two systems posess the same DSS critical solution. We therefore see 
that the Einstein-Maxwell theory is not the /i — * limit of the Einstein-Proca system, a discontinuity reminiscent 
of that observed in pure gravitation [14j . Indeed, since gravitation couples to all forms of energy, it couples to 
the longitudinal mode of the Proca field, amplifying it during spherically symmetric critical collapse relative to the 
transverse modes which become negligible. The physics of the critical gravitational collapse of a Proca field therefore 
becomes indistinguishable from that of a massless scalar. 

We close by noting that the type I critical solution we have found is essentially an analog for the Proca system of 
the soliton solution found by Seidel and SuenQ for the massive scalar field. We therefore expect that our solution 
could be found directly using the methods of reference^. We will address this issue in a separate paper. 
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